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Short  Term  Forecasting  of  Cloud  and  Precipitation 

Along  Communication  Paths 


I.  INTRODUCTION 

Communications  between  satellite  and  ground  station  systems  are  affected  by 
the  intervening  meteorological  environment.  With  the  trend  towards  higher  com¬ 
munication  frequencies,  knowledge  of  the  state  of  the  atmosphere  along  the  commu¬ 
nication  path  and  the  ability  to  provide  short  term  forecasts  of  this  environment, 
becomes  essential.  This  requirement  demands  knowledge  of  the  present  and  future 
three-dimensional  structure  of  precipitation  and  cloud  content  throughout  the  region 
of  interest.  This  report  is  concerned  with  a  review  of  data  sources  for  the  real¬ 
time  observation  of  cloud  and  precipitation  systems,  and  techniques  for  combining 
and  extrapolating  these  data,  to  prepare  the  desired  forecast.  There  are  two 
sources  of  data  which  provide  a  reasonable  data  base  for  solution  of  this  problem. 

First,  radar,  preferably  Doppler  Radar,  at  non -attenuating  frequencies  (for 
example,  3  GHz)  provides  the  best  mechanism  for  interrogating  the  interior  pre¬ 
cipitation  structure  of  storm  regions.  The  primary  radar  parameter  is  the  equiva¬ 
lent  radar  reflectivity  factor,  which  can  be  interpreted  as  a  measure  of  precipitation 
intensity  and  thus  can  be  related  to  attenuation  of  a  transmitted  signal.  Occasion¬ 
ally,  however,  the  radar  coverage  is  not  adequate.  Depending  upon  the  radar 
sensitivity,  a  certain  mass  of  precipitation  is  required  to  provide  sufficient  signal 
to  make  reliable  measurements.  Thus  the  radar  does  not  identify  the  total  storm 

(Received  for  publication  31  December  1985) 
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outer  boundary.  Also,  near  storm  top,  due  both  to  sensitivity  and  radar  scan 
requirements,  loss  of  information  may  also  occur.  Finally,  at  lowest  levels  where 
active  meteorological  events  such  as  gust  fronts  and  newly  initiated  convection 
may  be  present,  the  radar  observations  may  be  deficient.  Fortunately,  these 
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I  regions  where  radar  measurements  can  be  lacking  may  usually  be  successfully 

i  interrogated  by  satellite.  The  visible  and  infrared  imagery  data  are  highly  attenu¬ 

ated  in  storm  precipitation  regions,  making  the  satellite  well  suited  for  monitoring 
storm  tops,  the  outer  storm  environments,  and  regions  well  removed  from  active 
precipitation  generation.  The  radar  and  satellite  data  are  thus  highly  complementary 
and  may  provide  sufficient  and  appropriate  information. 

Fast  efforts  have  usually  been  interested  in  providing  forecasts  of  large  scale 
meteorological  features.  This  generally  resulted  in  use  of  large  time  (hrs)  and 
spatial  (tens  of  km)  scales,  or  the  use  of  a  single  parameter  to  represent  the  total 
storm  (for  example,  center  of  mass)  environment.  Here,  however,  forecasts 
are  desired  on  time  scales  of  5  to  30  min  and  spatial  scales  of  a  few  km.  Thus  this 
problem  is  significantly  different.  In  the  ensuing  discussions,  traditional  techniques 
for  storm  detection  and  tracking  will  be  presented.  They  represent  some  of  the 
basic  methodologies  still  in  use  today,  and  results  of  such  efforts  are  pertinent  to 
the  present  study.  Then  the  focus  will  be  directed  to  a  more  detailed  discussion  of 
the  current  problem,  and  proposed  methods  for  its  solution.  Last,  a  short  graphical 
demonstration  outlining  the  solution  methodology  will  be  given. 

2.  TRADITIONAL  TRACKING  TECHNIQUES 
2. 1  Radar  Data  Fields 

The  use  of  radar  data  for  storm  identification  and  tracking  has  been  a  subject 
of  interest  for  some  time.  The  feature  generally  monitored  has  been  the  radar 
reflectivity  factor.  The  general  manner  in  which  it  has  been  used  may  be  demon¬ 
strated  with  a  few  examples. 

Austin  and  Bellon  developed  the  Short-Term  Automated  Radar  Prediction 
(SHARP)  method  which  utilizes  simple  cross -correlation  analysis  of  the  entire  data 
field  to  determine  the  average  translation  of  the  precipitation  field  with  time.  This 
method  generally  performs  well  in  non -evolutionary  environments,  where  initiation 
and  dissipation  effects  are  slight,  where  the  environmental  winds  do  not  vary  greatly 
with  height,  and  where  the  precipitation  environments  are  shallow.  This  method 


1.  Bellon,  I).  R.  ,  and  Austin.  G.  L.  ( 1976)  SHARP  (Short-Term  Automated  Radar 
Prediction)  A  Real  Time  Test,  17th  Conference  on  Radar  Meteorology, 
Amer.  Meteorol.  Soc.  ,  Boston,  pp  522-525. 
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is  most  useful  when  no  well  defined  storm  structure  such  as  localized  cells,  or 
when  data  extends  beyond  the  viewing  region  of  analysis,  occurs.  This  method 
has  been  used  successfully  to  make  2 -hr  forecasts  of  large  scale  precipitation 
motions. 

2 

Elvander  compared  a  number  of  tracking  and  forecasting  techniques.  These 

included  the  SHARP  method,  a  linear  least  squares  (LLS)  tracking  of  individual 

centers  of  area  or  mass,  and  a  more  robust  technique  developed  by  the  Stanford 

3 

Research  Institute  (SRI)  (Wolf  et  al  )  which  employs  a  global  track  of  the  entire 
field  but  also  includes  adjustments  for  individual  storms  for  deviations  from  the 
field  average  value.  The  methods  were  employed  with  both  radar  reflectivity 
factor  and  vertically  integrated  liquid  (VIL)  data  as  input.  Whereas  the  SHARP 
method  results  in  one  translation  and  forecast  vector  for  the  entire  data  field,  the 
LLS  and  SRI  methods  return  track  and  forecast  vectors  for  the  individual  storms 
of  interest.  For  reflectivity  tracking  only  the  0  degree  elevation  data  were  used. 

In  most  situations,  including  large  and  small  distributions  of  storm  cells,  and 
frontal  and  squall  line  situations,  the  simple  cross -correlation  (SHARP)  method 
performed  best  in  forecasting  motions  of  individual  storms.  The  linear  least 
squares  method  was  least  successful.  In  tracking  VIL  the  LLS  method  was  pre¬ 
ferred.  Some  additional  observations  common  to  most  tracking  and  forecasting 
methods  were  found.  Tracking  and  forecasting  with  VIL  was  found  more  accurate 
than  with  use  of  reflectivity  factor  data  at  any  single  height  level.  The  track- 
ability  of  storm  cells  generally  increased  with  increasing  precipitation  intensity 
(reflectivity).  Forecast  errors  increased  with  forecast  lead  time  (here  a  linear 
trend  was  observed.  It  is  informative  to  interpret  these  results  for  some  apply 
generally  to  all  tracking  methodologies. 

The  relative  behavior  of  the  methods  when  using  VIL  vs  reflectivity  factor 
data  suggests  that  a  vertically  averaged  feature  is  easier  to  track  than  features  at  a 
single  height  level.  This  is  understandable  when  one  considers  that  at  a  single 
height  level  the  radar  contours  change  in  response  to  two  basic  effects,  evolution  and 
movement  of  precipitation  features  through  this  level  from  adjacent  levels.  This 
latter  effect  is  particularly  noticeable  when  a  coarse  horizontal  grid  (here  7  <  7  km) 
is  used.  As  the  precipitate  moves  up  and  down  within  the  storm  it  may  remain  at 
the  same  horizontal  grid  location  for  numerous  height  levels,  but  jumps  in  and  out 
of  these  same  levels  individually  as  the  precipitate  travels  vertically.  Thus  VIL. 
when  used  with  a  coarse  grid  is  relatively  insensitive  to  the  effects  of  vertical 

2.  Elvander,  R.C.  (107  6)  An  Evaluation  of  the  Relative  Performance  of  Three 
Weather  Radar  Echo  Forecasting  Techniques,  17th  Conference  on  Radar 
M  eteorologv,  Anter.  Meteorol.  Soc.  .  Boston,  pp  526-522. 

Wolf,  D.  F..  ,  Hall,  I).  .1.  ,  anti  Endlieh,  R.  M.  (  1077)  Experiments  in  automati 
cloud  tracking  using  SMS -GOES  data.  .1.  Appl.  Meteorol.,  16:1210-1230. 
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displacement  of  the  precipitate,  and  behaves  as  a  conservative  quantity.  The 
structure  at  a  single  height  level,  however,  may  vary  significantly.  In  tracking 
VIL  the  LLS  method  was  preferred  over  whole  field  cross -correlation  presumably 
because  the  individual  storm  cells  had  varying  propagation  velocities,  different 
from  the  average  field  value.  When  tracking  reflectivity  at  a  single  level  the  LLS 
method  was  less  accurate  than  the  SHARP  method  because  the  apparent  varying 
motions  induced  upon  the  center  of  mass  due  to  the  precipitate  vertical  motions 
and  vertical  wind  shear  were  greater  than  the  displacement  differences  between 
the  whole  storm  motions  and  the  field  average  value. 

The  SRI  method  determined  both  the  global  and  individual  storm  motions.  This 
method  pre-processed  the  input  data  to  describe  the  data  field  in  terms  of  clusters, 
where  clusters  consisted  of  neighboring  regions  of  similar  structure.  Cross¬ 
correlation  was  performed  on  the  clusters  to  determine  their  displacements.  As 
a  result  of  the  clustering  process,  this  method  was  most  successful  in  monitoring 
the  dissipation,  splitting,  and  merging  of  storm  cells.  However,  the  SRI  method 
was  computationally  intensive  and  offered  no  increase  in  accuracy  in  actual  tracking 
ability. 

A  slightly  modified  version  of  the  centroid  tracking  method  was  employed  by 
4 

Collier  where  comparisons  were  performed  between  whole  field  cross -correlation 
and  a  modified  storm  center  of  mass  tracking  technique.  Here  a  very  coarse  grid 
of  20  X  20  km  was  used.  The  data  were  pre-processed  using  a  low  threshold  to 
classify  areas  as  containing  rain  or  no  rain.  Adjacent  grids  of  rain  were  combined 
into  clusters.  The  centers  of  mass  of  the  clusters  were  determined  without  regard 
to  actual  precipitation  intensity,  thus  areal  centers  of  mass  were  determined.  These 
centers  were  then  individually  tracked.  Analysis  indicated  that  for  situations  where 
little  evolution  is  ongoing  and  the  echoes  are  moving  rather  uniformly,  whole  field 
cross -correlation  techniques  were  the  simplest  to  implement  and  quite  reliable.  In 
non-stationary  environments,  however,  tracking  the  individual  clusters  through  use 
of  the  center  of  mass  method  was  preferred.  Skilled  subjective  analysis  resulted 
in  very  little  increase  in  accuracy  of  the  forecasts  over  the  cluster -centroid  tech¬ 
nique.  Here  the  large  grid  filtered  out  evolutionary  and  vertical  transport  effects 
noted  in  the  previous  study.  This  allowed  for  individual  storm  tracking  to  be 
superior  over  applying  the  whole  average  field  motion  to  the  individual  storms  com¬ 
prising  the  field.  Thus  the  usefulness  of  these  various  tracking  algorithms  is  seen 
to  be  dependent  not  only  upon  the  meteorology,  but  also  upon  such  parameters  as 
grid  scale  size  in  the  horizontal  and  vertical  directions. 

4.  Collier,  C.G.  ( 198 1)  Objective  Rainfall  Forecasting  Using  Data  from  the 
United  Kingdom  Weather  Radar  Network,  Proc.  AamaP  Symposium, 


Further  assessment  of  short  term  forecasting  of  precipitation  using  radar 
data  was  considered  by  Ciccione  and  Pircher.  ®  Once  again  cross -correlation  and 
centroid  methods  were  studied,  but  with  a  horizontal  grid  scale  resolution  which 
could  be  reduced  to  2  km.  Both  methods  determined  the  whole  field  translation  with 
no  concern  for  the  deviating  motions  of  the  individual  storms  comprising  the  field. 

T  hey  noted  large  scale  phenomena  were  better  tracked  and  forecast  when  a  large 
grid  resolution  was  used.  Accuracy  decreased  with  increasing  forecasting  period 
and  the  standard  error  of  displacement  increased  from  zero  at  zero  lag  to  an 
asymptotic  value  near  20  percent  at  2  hrs  time  lag.  Success  also  increased 
with  increasing  proportion  of  rain  containing  area  and  with  increasing  verification 
area.  On  the  other  hand,  during  periods  where  the  rain  was  distributed  in  small 
areas,  cell  lifetimes  were  small,  and  grid  resolution  small,  success  increased 
with  decreasing  lag  times  and  small  resolution.  This  last  effect  simply  indicates 
that  as  the  pace  for  evolving  features  increases,  finer  spatial  and  temporal  resolu¬ 
tion  is  required.  Cross -correlation  was  clearly  the  superior  method.  The  results 
here  suggest  that  the  overall  motion  of  the  sum  of  the  cells  was  fairly  well  aligned 
with  the  field  motion,  but  that  a  fair  amount  of  deviation  of  individual  storms  from 
the  field  average  value  resulted  in  a  somewhat  more  erratic  track  for  the  field 
center  of  mass,  resulting  in  less  accuracy  when  field  center  of  mass  was  employed. 
This  result  once  again  underscores  the  difficulty  in  applying  whole  field  translation 
P  to  individual  storms. 

A  significant  ongoing  tracking  and  forecast  effort  is  that  being  employed  with 

6 

the  Next  Generation  Weather  Program  (NEXRAD).  Here  a  mass  weighted  centroid 
method  is  being  employed  for  prediction  of  future  individual  storm  locations,  as 
determined  by  the  storm  center  of  mass.  This  parameter  is  obtained  by  correlating 
the  various  layer  centers  of  mass  as  determined  from  data  acquired  from  a  sequence 
of  elevated  radar  scans  through  the  storm  of  interest.  Assuming  the  storm  centroid 
a  conservative  quantity  between  radar  scan  sequence  intervals,  the  displacement  of 
the  storm  center  of  mass  is  measured  through  use  of  a  nearest  neighbor  algorithm. 

A  linear  least  squares  method  is  applied  to  the  history  of  centroid  locations  to 
develop  forecasts  of  future  positions.  Although  a  useful  storm  tracking  and  fore¬ 
casting  technique,  it  offers  no  insight  into  the  evolving  nature  of  internal  storm 
structure. 


5.  t'iccone,  M.  ,  and  Pircher,  Y.  (1984)  Preliminary  Assessment  of  Yerv  Short 
Term  Forec  asting  of  Rain  from  Single  Radar  Data.  Proc.  Second  International 

•  Svmposium  on  Xoweasting,  Xorrkoping,  Sweden,  pp  24  1  -246. 
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This  short  summary  attempts  to  demonstrate  that  primary  efforts  in  precipita¬ 
tion  forecasting  have  been  focused  on  storm  scale  through  large  scale  phenomena. 
Large  scale  efforts  generally  track  only  gross  features  at  one  altitude,  while  storm 
scale  efforts  have  been  reduced  to  a  study  of  storm  centroid  motions.  These 
methods,  successful  in  their  own  applications,  do  not  address  the  current  problem. 
For  forecasting  cloud  and  precipitation  structure  along  a  communication  path,  it  is 
the  internal  three-dimensional  storm  structure  which  must  be  observed,  understood, 
and  forecast.  Thus  this  effort  requires  a  different  approach. 

2.2  Satellite  Data  Fields 

Analysis  of  satellite  data  has  essentially  followed  along  the  same  route  as  that 
taken  for  radar,  with  cross -correlation  and  centroid  methods  predominating.  Some 
potential  problems  associated  with  use  of  radar  data  were  briefly  noted.  Use  of 
satellite  visible  and  infrared  imagery  also  present  some  difficulties.  The  satellite 
frequencies  of  the  standard  GOES  type  satellite  are  generally  attenuated  by  inter¬ 
vening  cloud  and  precipitation.  As  a  result  of  this  only  the  outer  cloud  boundary 
that  faces  the  viewing  direction  of  the  satellite  is  observed.  Difficulty  in  differen¬ 
tiating  between  storm  side  and  top  may  result.  Use  of  complementary  infrared  data 
may  remove  some  of  this  ambiguity,  however  this  may  be  difficult  owing  to  the 
considerably  larger  footprint  of  the  infrared  image  as  compared  to  the  visible. 
Further  complications  arise  from  variation  in  observed  cloud  brightness  between 
observations.  This  is  partly  due  to  the  varying  orientation  of  the  storm  feature 
relative  to  the  satellite  and  the  sun,  both  of  which  may  change  differently  with  time. 
This  effect  may  again  result  in  false  identification  of  storm  features.  The  resulting 
outcome  of  these  concerns  is  an  ambiguity  in  the  vertical  and  horizontal  locations 
of  exterior  storm  features  relative  to  a  coordinate  system  fixed  to  the  earth's  sur¬ 
face.  Obviously  such  uncertainty  can  introduce  significant  error  in  the  placement 
of  tracked  features,  and  ultimately  significant  error  in  their  forecast  locations.  A 
further  difficulty  lies  in  the  initial  registration  of  the  entire  satellite  data  field 
itself.  Thus  with  use  of  satellite  data,  much  effort  is  taken  in  first  determining  the 
actual  location  of  features  of  interest,  both  in  height  and  horizontal  direction,  prior 
to  its  use  as  input  into  a  forecasting  technique. 

3 

As  an  example  of  a  very  robust  method  the  SRI  (Wolf  et  al  )  technique  discussed 
earlier  will  be  presented  in  some  greater  detail  as  it  has  been  used  with  satellite 
data.  The  first  task  is  landmarking.  This  is  accomplished  by  correlation  tech¬ 
niques.  First,  clouds  are  thresholded  out  bv  eliminating  brightness  values  above 
a  certain  threshold.  A  landmark  template  is  then  passed  over  the  imagery  field  in 
a  coarse  fashion  to  determine  rough  displacement  of  the  current  observation  from 
the  previous.  Next,  a  fine  template  is  employed  with  the  best  match  defining  the 


the  proper  field  location  in  geographical  sense.  Accounting  for  the  variation  in 
cloud  brightness  between  observations  is  done  statistically,  by  forcing  the  bright* 
ness  frequency  distributions  for  the  different  observational  periods  to  match  as 
well  as  possible.  The  alternative  approach  would  be  to  rely  upon  the  sinusoidal 
variation  of  sun  height  alone.  Once  properly  registered,  a  more  complete  dis¬ 
crimination  between  cloud  and  ground  features  is  performed  by  use  of  a  series  of 
thresholding  applications.  After  elimination  of  ground  targets  the  data  are  smoothed 
and  thresholded  again  to  retain  only  the  brighter  features. 

Clusters  are  now  formed  bv  estimating  the  brightness  deviations  between 
neighboring  grid  points.  Locations  where  sharp  deviations  are  detected  are  re¬ 
moved.  This  has  the  effect  of  locating  and  removing  the  outer  edge  of  the  actual 
cloud  boundaries,  white  still  retaining  the  shape  of  the  cloud  boundary.  Touching 
groups  are  now  formed  from  touching  regions  exhibiting  similar  brightness  values. 
Further  smoothing  and  elimination  of  small  isolated  features  to  be  considered  as 
noise  is  continued,  resulting  in  a  field  of  cloud  areas  which  form  the  tracers. 

Tracer  translation  is  done  either  by  pattern  recognition  or  correlation.  The  pattern 
recognition  method  compares  attributes  of  the  individual  tracers  (clusters)  to  deter¬ 
mine  the  translation.  Such  attributes  are  size,  location  of  center,  average  bright¬ 
ness,  average  IR  temperature,  and  the  rms  dispersion  of  the  tracer  region  along 
the  x  and  y  axes.  The  relative  translation  of  two  features  which  have  similar 
attribute  values  and  which  lies  within  predetermined  bounds,  such  as  mean  field 
translation  plus  a  deviation,  is  a  measure  of  the  displacement.  Cross -correlation 
of  successive  observations  of  individual  clouds  may  also  be  performed  to  determine 
the  different  translation  velocities.  In  any  case,  tests  are  performed  over  individual 
tracers  resulting  in  a  field  of  varying  tracer  motions.  Through  the  clustering 
process  this  method  provides  additional  information  concerning  storm  dissipation, 
splitting,  and  merging.  As  stated  earlier,  however,  it  does  not  necessarily  provide 
better  tracking  and  forecast  capability  over  the  simpler  correlation  and  center  of 
mass  techniques  for  the  storm  and  field  scales.  Also,  as  currently  employed,  it 
still  does  not  determine  the  evolution  of  distinct  structures  such  as  contour 
boundaries. 

The  methods  discussed  may,  varying  with  meteorological  situation,  return  a 
reasonable  estimate  of  cloud  or  precipitation  feature  motion  with  time  if  this 
feature  is  not  evolving  on  a  time  scale  comparable  to  a  few  observation  periods. 

The  problem  being  addressed  here  demands  more  attention  to  the  evolutionary  nature 
of  the  observed  systems  as  well  as  the  general  translation  of  the  features.  This 
can  be  accomplished  only  by  monitoring  the  motion,  shape,  size,  and  intensity  of 
cloud  and  precipitation  structures  with  time.  Of  the  few  techniques  discussed  here, 
the  SHI  method  is  most  closely  aligned  to  the  problem  as  it  is  seen  here.  It  offers 


the  potential  to  provide  a  measure  of  the  change  in  shape  and  size  required  to 
accurately  forecast  the  trend  in  the  tracked  feature.  However,  the  SRI  model 
removes  significant  amounts  of  actual  data  and  is  a  heavy  computational  burden. 

A  number  of  other  methods  exist  for  monitoring  the  spatial  distribution  of  the  cloud 
and  precipitation  data.  Spectral  techniques  may  be  employed  to  analyze  the  distribu¬ 
tion.  However,  the  spectral  information  does  not  easily  lend  itself  to  determining 
the  shape  and  size  at  any  single  location.  On  the  other  hand,  shape  fitting  methods 
and  boundary  tracking  techniques  appear  to  be  useful  alternative  methods  for 
monitoring  the  changes  of  shape  and  size  of  the  tracers.  The  remainder  of  the 
discussion  will  focus  on  their  application  to  the  present  problem. 

3.  TIIE  FORECAST  PROBLEM 

The  goal  ofihis  effort  is  to  provide  short  term  forecasts,  for  5-  to  10-min 
intervals  for  a  maximum  length  of  30  min,  with  a  spatial  resolution  of  a  few 
km.  This  can  be  accomplished  only  by  monitoring  the  motion,  shape,  size,  and 
intensity  of  cloud  and  precipitation  structures  with  time.  A  compromise  must  be 
struck  between  the  spatial  scales  and  the  grid  scale  employed  to  accommodate 
the  expected  5-  to  10-min  update  interval  of  radar  data  and  potential  storm 
feature  evolution.  To  forecast  for  a  given  location  we  must  have  some  knowledge 
>f  the  precipitation  content  within  a  region  which  could  influence  that  location  by 
simple  advection  or  other  means  over  the  maximum  predictive  time  interval.  Thus 
we  need  to  know  the  translation  speed  and  direction  for  the  data  in  an  appropriate 
corridor.  If  this  location  is  to  be  truly  arbitrary,  then  we  need  to  be  able  to  pro¬ 
ject  the  vector  displacement  for  the  entire  field  of  data.  Simple  field  translation 
may  be  sufficient  for  very  short  forecast  periods,  over  which  the  field  is  nearly 
stationary  (5  to  10  min)  in  structure.  However,  for  longer  forecasting  intervals 
(20  to  30  min),  the  structure  probably  will  not  be  effectively  stationary.  Evolution 
at  a  given  height  level,  whether  a  result  of  real  precipitate  distribution  growth  at  that 
level,  or  a  response  of  passage  through  this  level  from  an  adjacent  level,  should  be 
expected. 

Proper  solution  of  this  problem  requires  knowledge  of  the  behavior  of  the  pre¬ 
cipitation  structures  at  various  height  levels  in  the  environment.  It  was  noted  that 
tracking  vertical  integrated  liquid  resulted  in  smaller  errors  than  when  tracking 
features  at  a  single  height  level.  Thus  some  form  of  vertical  averaging  is  desired. 
However,  for  projection  of  cloud  and  precipitation  amount  along  an  arbitrary 
communication  path,  some  knowledge  or  the  vertical  structure  is  required.  This 
is  easily  seen  in  storm  regions,  such  as  with  the  anvil  or  high  reflectivity  factor 
aloft,  where  the  precipitation  content  there  should  not  be  attributed  to  all  height 
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levels  as  may  occur  with  a  vertically  averaged  quantity.  Another  reason  to 
retain  some  vertical  structure  is  that  knowledge  of  vertical  motion  of  a  precipita¬ 
tion  volume  provides  some  forewarning  for  new  development  at  adjacent  height 
levels.  Thus  new  development  at  a  lower  (or  higher)  level  need  not  always  be 
unexpected,  and  a  feedback  mechanism  is  available.  Thus  for  a  proper  identifica¬ 
tion  of  cloud  and  precipitation  along  some  arbitrary  path,  some  facility  for  retention 
of  vertical  structure  is  desired. 

In  this  effort,  description  of  the  field  of  precipitation  or  cloud  intensity  will  be 
accomplished  by  use  of  radar  reflectivity  factor  for  radar  data  and  brightness  level 
for  satellite  visible  imagery.  These  parameters  are  most  easily  transformed  into 
measures  of  attenuation.  The  data  field  structure  at  any  single  height  level  will 
consist  of  a  series  of  intensity  levels  and  will  appear  in  the  field  image  as  a  series 
of  contour  boundaries,  some  closed,  some  open.  These  contours  will  be  termed 
the  "features"  of  the  meteorological  field.  These  features  may  be  described  in  a 
number  of  ways  through  use  of  descriptive  parameters  which  will  be  termed 
"attributes".  Thus  in  the  final  analysis,  tracking  and  forecasting  of  a  field  of  data 
will  ultimately  result  in  tracking  or  forecasting  a  set  of  attributes.  A  demonstration 
of  the  use  of  attribute  parameters  in  the  forecasting  problem  is  presented  in 
Section  7. 

4.  FEATURE  DEFINITION  AND  EXTRACTION 
4.1  Field  Analyaia  vs  Comparison 

The  basic  methods  for  feature  extraction  are  generally  derived  from  (a)  field 
comparison  or  (b)  field  analysis.  Field  comparison  implies  taking  successive 
field  patterns,  or  subsets  of  field  patterns  at  different  times,  and  iteratively 
adjusting  the  location  of  one  field  relative  to  the  other  until  a  measure  of  maximum 
comparison,  usually  a  maximum  cross -cor relation,  is  obtained.  Here  the  feature 
of  interest  would  be  an  intensity  contour  boundary.  Quite  often  one  assumes  the 
structure  of  the  meteorological  field  is  stationary  and  evolution  is  of  no  concern. 

In  this  case  one  generally  employs  a  simple  form  of  field  comparison  to  determine 
the  simple  field  translation. 

Field  analysis  requires  some  analysis  of  the  data  and  its  subsequent  descrip¬ 
tion  in  terms  of  a  limited  number  of  descriptive  attributes.  When  the  new  field  is 
described,  then  a  comparison  is  made  between  the  two  sets  of  attributes.  Des¬ 
cription  of  the  fields  mav  include  such  parameters  as  centroid  locations,  where 
each  location  mav  represent  a  simple  areal  average,  or  an  intensity  weighted  value 
for  each  intensity  class.  Alternately,  the  descriptors  may  be  derived  from  some 


form  of  pattern  analysis.  In  any  case  the  primary  difference  is  that  field  analysis 
attempts  to  determine  the  underlying  patterns  of  the  data  field  whereas  comparison 
develops  no  understanding  of  these  basic  underlying  structures. 

Pattern  analysis  may  be  used  to  provide  a  description  of  the  contour  of  a  given 
intensity  class.  The  methods  may  include  fitting  the  contour  in  small  segments 
using  straight  lines  or  arcs.  Alternatively  the  actual  shape  of  the  entire  contour 
may  be  fit  by  polynomial  or  polygonal  means,  or  by  combinations  of  basic  shapes 
such  as  ellipsoids  as  utilized  in  vortex  analysis.  Such  methods  offer  different 
degrees  of  complexity  for  generation  of  the  attributes,  and  in  their  use  in  deter¬ 
mining  the  contour  changes.  In  the  forecasting  effort  here,  detection  of  field 
changes  and  their  forecast  effects  are  required.  Thus  here  we  emphasize  field 
analysis.  The  attributes  describing  these  features  may  vary  from  codes  defining 
line  segments  to  the  constants  forming  the  polynomial  fit  to  the  contour  boundary. 

Description  of  the  field  changes  may  be  accomplished  in  a  number  of  ways 
once  the  general  translation  between  two  successive  intervals  is  determined.  In 
field  analysis,  one  may  forecast  by  noting  the  changes  in  the  field  attributes.  A 
forecast  may  be  derived  by  determining  the  evolution  of  the  attributes  and  fore¬ 
casting  their  future  values.  Alternatively,  one  may  forecast  the  change  in  the 

attributes  and  adds  these  to  the  initial  attribute  states.  In  vortex  analysis,  as 
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demonstrated  by  Williamson,  a  contour  is  fit  by  the  resulting  boundary  of  over¬ 
lapping  ellipsoids  of  varying  sizes  and  orientation.  Thus  this  method  may  include 
forecasting  of  basic  ellipsoid  parameters  such  as  amplitude  and  center.  Shape 
fitting  requires  a  wholly  contained  boundary.  Contour  fitting  may  be  accomplished 
on  a  partial  boundary.  Thus  the  applicability  of  pattern  analysis  may  vary  with 
field  structure. 

The  current  effort  requires  detection  of  field  changes  and  development  of  their 
time  histories  for  forecast  purposes.  Particular  methods  for  feature  extraction 
will  be  discussed.  One  of  the  major  constraints  of  this  effort  is  the  rapid  data  and 
forecast  update  rate  of  5  to  10  min.  This  places  severe  limitations  on  the 
types  of  analysis  that  may  be  performed,  and  the  definition  of  attributes  by  which 
the  features  of  the  field  will  be  tracked.  Noting  that  the  basic  image  size  involved 
in  this  analysis  is  512  by  512  pixels,  and  that  within  each  image  there  will  likely  be 
a  number  of  features  to  be  monitored  it  seems  highly  unlikely  that  analyses  can  be 
performed  pixel  by  pixel  on  a  routine  basis.  In  addition  a  number  of  images  may  be 
required  to  be  retained,  thus  making  storage  of  field  data  a  problem.  Focus  will 
now  be  directed  towards  developing  appropriate  and  efficient  methods  for  the 
description  of  field  attributes. 


7.  Williamson.  I).  I..  (1081)  Storm  track  representation  and  verification, 
Tellus.  38:518-530. 


4i  Feature  Detection  aid  Definition 


Feature  detection  will  be  a  result  of  some  form  of  picture  segmentation,  where 
this  is  broadly  defined  as  "picking  out  the  objects  of  interest  from  the  background". 
The  two  general  approaches  for  obtaining  the  features  are  structural  and  statistical. 
Structural  methods  are  to  be  employed  when  the  data  field  contains  features  that 
have  obvious  coherency,  or  patterns  discernible  by  eye.  This  is  the  expected  state 
of  affairs  in  typical  radar  data  as  shown  in  Figure  1,  where  a  complex  set  of 
quantized  reflectivity  factor  contours  are  evident. 

There  are  a  number  of  ways  of  detecting  these  boundary  contours.  The  pro¬ 
cedure  may  be  considered  three-fold,  namely;  (a)  silhouette,  where  the  data  field 
is  thresholded  with  all  pixel  locations  having  a  value  above  the  threshold  given  a 
prescribed  value  (for  example,  one),  and  all  pixels  with  values  below  the  threshold 
a  value  of  zero);  (b)  search,  where  the  first  boundary  point  of  the  contour  of 
interest  is  located;  and  (c)  follow,  where  the  contour  is  followed  and  described  in 
some  manner.  Assuming  structural  methods  are  in  use,  a  number  of  boundary 
detection  methods  are  available.  The  simplest  is  the  gradient  operator,  that 
estimates  the  data  gradient  from  the  pixel  locations  about  the  pixel  of  interest. 
Because  the  data  field  has  been  described  in  binary  fashion,  all  possible  gradient 
values  are  known  in  advance,  and  boundary  crossings  are  easily  determined.  More 
complex  transformations  exist  for  determining  boundaries,  but  a  simple  gradient 
operator  is  both  quick  and  simple  and  should  be  adequate. 

When  the  data  field  presents  a  noisy  appearance,  containing  no  easily  de- 
scribable  features,  statistical  methods  may  first  be  employed  to  locate  those 
regions  that  have  similar  structure  in  the  statistical  sense.  Here  the  field  of  data 
is  described  in  non-descriptive  terms,  such  as  the  form  of  the  intensity  distribution, 
its  mean  and  variance,  either  for  the  entire  field,  or  more  likely  for  subset  regions 
of  the  large  scale  field.  After  thresholding,  the  data  field  is  first  segmented  if 
there  are  filled  and  closed  areas.  Each  area  is  now  investigated  to  determine  its 
intensity  distribution  and  associated  parameters.  When  areas  are  determined  to 
have  differing  parameters,  a  new  feature  area  has  been  located.  Here  then  a 
contour  represents  a  boundary  between  two  regions  having  differing  distribution 
parameters.  This  edge  mav  then  be  traced  to  define  its  boundary  as  described  in 
the  pure  structural  methodology. 
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Figure  1.  Contours  of  Radar  Reflectivity  Factor 
at  3.  35  km  Altitude  for  Storm  Complex  Observed 
Near  Wallops  Island,  Virginia.  Contours  are  in 
recycling  dot,  dash,  solid  pattern.  Minimum  and 
maximum  contour  plots  shown  are  15  and  40  dBz, 
respectively.  Positions  are  relative  to  radar 
location 


4.3  Boundary  Description 

The  description  of  the  boundary  is  now  pursued.  It  is  assumed  that  the  bound 
ary  has  been  determined.  A  number  of  methods  for  description  are  available, 
from  very  simple  pattern  methods  to  complex  curve  fitting.  The  most  primitive 
methods  employ  icons  or  position  chain  codes.  Icons  may  be  fixed  length  vectors 
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in  the  eight  possible  directions  of  movement  from  the  current  image  location 
(Freeman  chain  code),  or  combinations  of  straight  lines  or  arcs.  An  example  of 
the  simplest,  albeit  a  very  attractive  technique,  is  the  Freeman  chain  code  icon 
representation  shown  in  Figure  2.  The  boundary  is  described  by  a  series  of  code 
characters  (for  example,  here  values  1  through  8).  Alternatively,  a  position  chain 
code  could  be  a  series  of  successive  (x,  y)  pixel  locations.  The  next  level  of 
complexity  would  be  fitting  boundary  segments  of  varying  length  by  straight  lines, 
arcs,  or  simple  polynomials.  Polygonal  or  polynomial  approximation,  or  vortex 
fitting,  to  the  full  boundary  would  represent  about  the  most  complex  form  of 
boundary  description  method. 


(a)  (b)  (c) 

CHAIN  CODE  6  7  77  I  781  I  I  3  3  3  5  43  3  5  5 

(d) 

Figure  2.  Description  of  Contour  Boundary  Using  Chain  Coding  With 
(a)  the  Eight  Freeman  Chain  Code  Elements,  (b)  the  Contour  to  be 
Described,  (c)  Following  the  Contour  With  Chain  Code  Elements, 
and  (d)  the  Resultant  Chain  Code  Stored  in  Memory 

There  are  inherent  tradeoffs  with  each  level  of  complexity.  The  primitive 
techniques  are  very  fast  computationally.  For  example,  the  Freeman  chain  code, 
if  used  judiciously,  and  employing  a  memory  of  the  previous  value,  results  in  less 
than  four  searches,  on  the  average,  to  find  the  next  chain  component.  Their  use¬ 
fulness  is  determined  by  the  difficulty  with  which  they  may  be  transformed  into  a 
form  suitable  for  higher  level  analysis.  For  example,  how  difficult  is  it  to  trans¬ 
form  the  code  back  into  a  spatial  representation  to  determine  differential  boundary 


change?  Advancing  to  higher  level  methods  the  computational  intensity  increases 
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nonlinearly  with  increasing  complexity  of  the  technique  employed.  Thus  polynomial 
curve  fitting  is  quite  computationally  intensive  when  compared  to  the  more  primitive 
methods.  However,  it  may  be  more  easily  transformed  into  the  spatial  form  and 
more  readily  useable  for  subsequent  analysis.  There  exists  a  concern  in  the  use 
of  high  level  boundary  description,  however.  For  accurate  forecasting  on  small 
spatial  scales,  as  much  information  about  the  contour  boundary  structure  as 
possible  must  be  retained.  Thus  notches,  cusps,  particular  features  of  the  bound¬ 
ary  which  make  it  unique  from  other  contours  in  the  data  field  must  be  retained. 

Shape  fitting,  through  use  of  high  level  routines,  invariably  produces  a  highly 
smoothed  contour  shape.  The  boundary  features  that  are  lost  are  precisely  those 
which  may  be  very  important  to  be  retained.  The  primitive  methods,  however, 
offer  potential  for  memorizing  all  features  in  the  boundary.  For  example,  the 
Freeman  chain  code  logs  the  contour  pixel  to  next  pixel,  thus  allowing  for  no  informa¬ 
tion  loss.  These  concerns  make  the  primitive  techniques  the  more  attractive, 
although  later  more  extensive  analysis  may  be  more  involved. 

Finally,  in  all  methods  some  facility  must  be  present  to  allow  for  data  field 
smoothing  and  data  rejection.  A  casual  observation  of  a  typical  radar  reflectivity 
factor  image  (Figure  1)  makes  the  most  casual  observer  aware  of  the  generally 
very  complex  intensity  structure.  Smoothing  operations  will  be  generally  per¬ 
formed  prior  to  the  feature  extraction  process  to  reduce  the  noisy  appearance, 
which  is  a  result  of  true  meteorological  complexity,  measurement  methods  and 
errors.  If  a  primitive  boundary  description  method  is  used,  it  is  much  more 
difficult  to  remove  noise  from  the  code  sequence  than  first  smooth  and  then  code  the 
contour.  If  a  high  level  method  such  as  polynomial  fitting  is  employed,  the  noise 
may  introduce  significant  errors  into  the  fit  relation.  Rejection  of  data  may  be 
done  by  the  requirement  for  a  minimum  areal  extent  enclosed  by  a  boundary  contour. 
This  assumes  that  noise  is  random  and  adjacent  pixel  values  will  not  uniformly  rise 
above  the  predetermined  threshold  for  a  contiguous  set  of  pixels.  The  requirement 
for  vertical  structure,  however,  does  force  interrogation  of  adjacent  levels  to  deter¬ 
mine  whether  a  region  is  part  of  a  precipitation  or  cloud  volume  extending  to  other 
levels,  in  which  case  areas  which  fall  below  the  areal  threshold  may  be  retained. 

5.  FEATURE  MOTION  AND  EVOLUTION 
5.1  Feature  Translation  and  Rotation 

Once  the  feature  boundaries  are  determined,  overall  motion  of  these  features 
must  be  determined.  This  may  be  performed  by  cross-correlation,  or  by  image 
processing  methods,  of  the  contour  boundaries.  In  cross-correlation  the  contours 


for  two  successive  observations  are  translated  and  rotated  until  the  minimum  mean 
square  error  is  obtained.  Here  error  is  determined  by  measuring  the  difference 
in  contour  boundary  locations  for  a  single  image  row,  and  then  summing  up  the 
squared  errors  for  all  rows.  Iteratively  moving  one  image  around  the  general 
location  of  the  other  until  this  value  is  a  minimum  yields  the  optimum  contour 
translation  and  rotation  components  and  minimizes  the  overlap  difference  between 
the  two  contours.  An  alternative  approach  to  find  the  condition  of  maximum  over¬ 
lap  is  to  perform  video  subtraction  of  the  two  images  after  each  iteration  until  the 
minimum  area  results.  The  results  of  these  two  different  processes  are  the  same. 
Image  video  subtraction  and  area  sum  can  be  performed  very  efficiently  in  the 
image  processor,  and  is  the  preferred  route.  Determining  feature  motion  by 
measuring  the  displacement  of  the  center  of  mass  is  not  considered  appropriate 
since  it  does  not  allow  for  the  effects  of  rotation  of  the  boundary. 

S.2  Feature  Evolution 

Once  the  general  overall  motion  of  a  contour  has  been  obtained  the  change  in 
contour  structure  must  be  determined.  With  a  history  of  the  translation,  rotation, 
and  change,  the  future  feature  structure  may  be  forecast.  There  are  numerous 
methods  for  detecting  these  changes.  However  only  three,  which  demonstrate  three 
different  levels  of  complexity,  will  be  discussed. 

First,  assume  a  very  primitive  contour  descriptor  has  been  employed.  The 
Freeman  chain  code  will  serve  as  our  example.  Our  information  consists  of  the 
starting  pixel  location  for  the  contour,  and  a  series  of  code  characters.  Other 
information  such  as  overall  translation  and  rotation  are  also  assumed  available. 

One  may  compare  the  two  successive  chain  codes  to  locate  positions  along  the  chain 
where  the  codes  differ.  This  allows  one  to  immediately  focus  upon  those  areas 
where  measurable  change  in  boundary  structure  is  occurring.  Chain  segments 
where  the  codes  have  remained  the  same  are  assumed  to  have  no  detectable  change 
during  this  time  interval.  The  chain  code  sequences  of  interest  are  transformed 
into  spatial  terms  for  interpretation.  Alternatively  a  chain  code  grammar  has  been 
developed  to  eliminate  the  spatial  transformation  step  and  provide  for  immediate 
interpretation.  In  any  case,  with  the  changes  interpreted  in  a  spatial  context,  time 
histories  are  developed,  understood,  and  forecasts  of  future  change  determined. 

The  next  level  of  effort  employs  what  may  be  termed  a  raster  vector  approach, 
potentially  most  useful  for  closed  contour  regions.  W  ith  the  contour  data  residing 
in  the  image  processor  one  may  consider  defining  it  bv  listing  the  starting  and  end 
points  for  each  successive  horizontal  raster  line.  Thus  the  contour  will  be  a  chain 
code,  with  the  elements  of  the  code  being  the  element  line  number,  and  starting 
and  ending  points  (n,  xl,  x2>,  as  many  sets  as  required  to  describe  the  contour 


positions  along  that  single  raster  line.  In  the  initial  description,  all  raster  vectors 
will  be  horizontally  oriented.  A  succession  of  observations  of  the  contour  of 
interest  will  result  in  a  series  of  measurements  of  its  translation,  rotation,  and 
chain  code  sequences.  One  may  choose  to  follow  each  raster  vector  element  in 
time.  For  this  case  the  chain  code  elements  must  be  described  on  a  grid  which  is 
allowed  to  translate  and  rotate,  thus  applying  the  translation  and  rotation  to  each 
raster  vector.  In  this  manner  one  is  able  to  measure  the  evolution  of  a  particular 
set  of  contour  relative  locations,  either  as  pair  or  individually,  to  develop  their 
histories,  and  forecast  their  future  change.  Alternatively,  one  may  define  raster 
vectors  as  always  being  horizontally  oriented,  even  though  the  contour  may  be 
rotating  in  time.  In  this  case  a  time  history  of  a  raster  vector  element,  a  given 
raster  vector  number  (n),  may  eventually  represent  different  locations  along  the 
contour.  In  this  method,  which  will  be  termed  the  "simple"  raster  vector  method, 
the  rotation  is  actually  accounted  for  by  the  raster  vector  end  point  variations. 

Both  methods  will  be  tested  to  determine  their  applicability  to  this  problem. 

As  a  final  example,  it  is  assumed  that  a  shape  fitting  technique  has  been 
employed  to  describe  a  closed  contour.  The  vortex  fitting  routine  will  serve  as 
our  example.  A  contour  has  been  described  by  a  series  of  overlapping  ellipsoids 
of  varying  amplitudes,  centers,  eccentricities,  and  other  attributes.  With 
successive  observations  a  sequence  of  ellipsoid  attributes  is  developed.  With  a 
time  history  of  attributes  for  each  individual  ellipsoid,  future  ellipsoid  shapes, 
sizes,  orientations,  and  locations  may  be  forecast.  This  composite  set  of  forecast 
ellipsoids  then  represents  the  forecast  contour  in  location,  shape,  and  size.  In 
this  approach  actual  measurement  of  feature  translation  and  rotation  is  not  neces¬ 
sary,  since  it  is  essentially  being  accounted  for  by  variation  in  the  attributes  of 
the  ellipsoids  describing  the  contour.  Feature  translation  and  rotation  may  be 
determined  aposteriori  if  desired. 

It  has  been  shown  how  the  variables  describing  a  contour  of  interest  may  be 
tracked  and  their  time  histories  developed  for  use  in  making  forecasts.  Alterna¬ 
tively,  instead  of  tracking  the  time  varying  attributes,  the  change  in  the  attributes 
may  be  tracked,  and  forecasts  built  upon  adding  to  an  initial  reference  attribute 
value  a  forecast  of  the  change  of  the  attribute  from  this  reference,  l  or  example, 
in  the  vortex  fitting  scheme,  instead  of  tracking  ellipsoid  amplitude,  we  make  the 
initial  amplitude  the  reference  value  and  track  the  change  in  the  amplitude  from  the 
reference  value.  This  method  has  the  advantage  of  allowing  one  to  use  the  second 
time  derivative  of  the  quantity  of  interest  in  the  forecast  process,  since  the  quantity 
being  tracked  is  the  change  in  the  attribute  rathei  than  the  attribute  itself.  W  hether 
this  is  truly  beneficial  depends  somewhat  upon  the  meteorology  of  the  situation,  the 
rate  of  feature  evolution,  the  noise  in  the  data,  and  the  method  used  for 
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feature  description.  Another  factor  which  must  be  considered  is  that  tracking 
change  in  a  variable  rather  than  the  variable  itself  results  in  less  information. 

Two  successive  locations  yield  two  attribute  values  but  only  one  attribute  change 
value.  It  is  quite  possible  that  the  lifetimes  of  the  features  may  be  so  short  as  to 
demand  use  of  as  much  data  as  possible  to  remove  the  effects  of  noise  to  discern  the 
true  trend  and  thus  make  a  reasonable  forecast.  This  would  preclude  use  of  the 
change  method.  In  the  final  analysis  both  methods  will  be  tested  to  determine  which 
allows  for  the  most  accurate  forecasts. 


6.  FORECAST  TOOLS 

The  method  by  which  a  forecast  is  made  is  somewhat  dependent  upon  the 
feature  description  approach  employed.  If  a  relatively  primitive  descr  ption  tool 
such  as  chain  coding  is  used,  then  boundary  description  processing  is  very  quick. 
Here,  however,  a  large  number  of  elements  may  be  required  to  develop  the  full 
contour  code,  and  to  forecast  values  for  each  item  in  the  data  set  a  large  number 
of  computations  may  be  involved.  Thus  a  simple  forecasting  method,  requiring 
very  few  computations  and  little  computer  memory  allocation  is  desired  for  the 
primitive  methods.  Conversely,  a  high  level  contour  description  method  such  as 
vortex  fitting  requires  only  a  small  number  of  parameters  to  be  retained  and 
tracked.  The  potential  for  use  of  a  complex  forecasting  technique  is  now  present. 

It  must  be  noted,  however,  that  this  is  not  necessarily  a  clear  cut  possibility,  for 
the  time  involved  in  developing  the  contour  description  in  the  high  level  method  may 
preclude  use  of  any  but  the  simplest  forecasting  tools.  Only  experience  with  real 
data  analysis  will  provide  a  guideline.  In  this  effort,  forecasting  techniques  ex¬ 
hibiting  a  broad  range  of  complexity  will  be  tested.  However,  because  of  the  hard 
time  limit  available  for  data  and  forecast  updates,  it  is  expected  that  the  main  focus 
will  initially  remain  with  the  simpler  methods. 

For  the  moment,  concern  will  focus  only  on  simple  techniques,  and  little  con¬ 
sideration  will  be  paid  to  whether  the  variables,  or  their  changes,  are  being  tracked. 
Potential  applications  will  be  discussed  later.  The  simplest  forecast  model  of  any 
practical  use  is  the  simple  smoothing  filter.  This  filter  inherently  assumes  the 

data  to  be  stationary  since  it  does  not  include  a  trend  factor.  A  useful  form  is  the 

8 

simple  exponential  filter,  as  demonstrated  by  Makridakis  and  Wheelright,  whicb 
may  be  written  as 


8.  Makridakis,  S. ,  and  Wheelwright,  S.  C.  (19781  Forecasting  Methods  and 
Applications,  John  Wiley  and  Sons,  New  York. 
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(1) 


Ffl“Kxt+(1-K>rt  • 
where  K  is  a  weight  in  the  range  of  0  to  1.  This  relation  may  also  be  written  as 
l't+1  =  KX  f  K(l-K>Xt  l  . . .  K(l-K>N'IXt  N+1  (2) 

*tVl-r.*KEt  ,3> 

whe  re 


VXt' 


(4) 


In  these  relations  Xt  is  the  measured  variable  at  the  current  observation,  F{ 
is  the  previous  forecast  for  the  current  time,  E^  is  the  measured  error  of  the  fore¬ 
cast  for  the  current  observational  time,  and  F.  ,  is  the  forecast  for  one  time 

t+1 

interval  into  the  future.  Depending  upon  the  form  of  the  relation  one  may  interpret 
this  filter  [Eq.  (3)]  as  forming  a  new  forecast  by  adding  to  the  current  observation 
the  error  of  the  forecast  for  this  current  observation.  Equation  (2)  indicates  that 
this  filter  naturally  gives  greater  weight  to  current  observations  than  to  previous 
ones. 

Another  version  of  a  smoothing  filter  is  the  moving  average  filter,  which  is 
written  as 


r*i-(Vxt-i-  xtW/N 


(5) 


=  Xt/N-Xt-N+l/N+Ft 


(6) 


where  all  observations  are  given  the  same  weight.  As  noted  these  two  filters  may 
be  useful  relations  if  the  parameter  is  basically  stationary  and  variations  in  the 
parameter  reflect,  noise  input.  In  essence  these  filters  attempt  to  determine  the 
true  mean  value  by  smoothing  out  the  noise. 

At  first  glance  these  simple  routines  may  seem  of  little  use,  however,  if  we 
are  tracking  the  change  in  a  parameter  such  as  position,  then  we  are  essentially 
estimating  an  average  velocity  of  translation,  which  is  now  a  very  useful  parameter 
in  forecasting  quantities  such  as  motion  of  large  scale,  or  stratiform  environments. 
The  exponential  filter  employs  the  concept  of  negative  feedback  |Eq.  <3>),  where 


j 

I 
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the  new  forecast  is  adjusted  for  the  error  oncurred  in  the  previous  forecast.  Also, 

■  the  exponential  filter  allows  the  option  for  varying  the  relative  weight  given  to 

J  recent  vs  past  observations.  This  is  particularly  relevant  for  meteorological 

!  systems  where  feature  evolution  may  demand  that  the  history  of  the  parameter 

focus  primarily  on  recent  observations.  In  addition  the  exponential  filter  requires 
very  little  storage  memory,  with  only  the  present  observation  and  the  previous 
forecast  of  the  parameter  required  to  make  a  new  forecast.  The  moving  average 
filter,  on  the  other  hand,  requires  retention  of  N  observations  and  is  more  compu¬ 
tationally  intensive.  Thus  the  simple  exponential  filter  may  be  useful  when  used 
with  the  primitive  feature  description  methods  where  many  parameters  are  updated. 

As  stated  above,  however,  these  simple  smoothing  filters  inherently  assume 
the  variable  to  be  stationary.  For  the  exponential  form  one  also  must  decide  upon 
a  value  for  a  weighting  factor.  This  last  concern  may  be  removed  by  use  of  the 
adaptive  simple  exponential  filter  which  automatically  adjusts  the  weight  value 
dependent  upon  the  error  detected  in  the  last  forecast.  With  these  simple  filters 
forecasting  for  more  than  one  time  interval  ahead  requires  the  assumption  that 


Xt+2=Ft+2 


(7) 


Xt+3  Ft+3 


(8) 


and  so  on.  Thus  one  sees  that  forecast  values  tend  to  a  constant  value,  since  all 
forecasts  are  assumed  to  have  zero  error.  Again,  whether  or  not  this  is  a  serious 
drawback  may  be  dependent  upon  the  variable  being  forecast. 

The  next  obvious  step  is  to  consider  a  parameter  as  exhibiting  a  linear  trend. 
The  standard  linear  regression  technique  uses  a  formulism  of 


F,  =  A  *  BX, 
t+  m  t+  m 


(9) 


where  ^t+rn  is  the  independent  variable  (generally  time)  at  the  future  time  m  inter¬ 
vals  ahead,  and  F,  is  the  forecast.  The  parameters  A  and  B  are  determined 
t+m 

through  the  relations 

A  =  SF/N  -  BSX/N  (10) 

B  =  <N£FX-XF£X>/(\XX2  -  ClX  )2>  . 


(11) 
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It  is  noted  that  again  N  values  of  both  F  and  X  are  retained  in  memory  for  each  new 
computation,  and  computational  complexity  has  increased  significantly  over  the 
simpler  filters  previously  discussed. 

An  easier  way  to  incorporate  trend  in  the  forecast  process  is  by  use  of  a 
modified  form  of  the  exponential  or  moving  average  filter.  The  Holt  two-parameter 
linear  exponential  filter  is  of  the  form 


F 


t+m 


St  +  m 


(12) 


where 

Bt  =  A(St  -  St_j)  +  (1  -  A)  Bt  l  (13) 

S  =CX  +(1-C)(S  H+BH).  (14) 

In  this  formulation  Fq.  (14)  corrects  St  for  the  error  in  the  last  estimate  by 
adding  in  the  previous  trend  estimation,  while  Fq.  (13)  effectively  updates  the  trend. 
Thus  both  smoothing  of  the  data  and  linear  trending  is  accounted  for  in  this  method. 
A  potential  drawback  is  the  use  of  the  two  weighting  constants  A  and  C.  These 
terms  should  generally  be  determined  through  some  form  of  error  analysis  (for 
example,  mean  square  error  minimization).  This  requirement  could  seriously 
degrade  its  real  time  applicability.  If  a  useful  set  of  weights  applicable  to  many 
meteorological  situations  can  be  determined  through  analysis  of  case  studies,  then 
this  method  is  preferred  over  the  standard  linear  regression  formulism  for  the 
following  reasons.  The  method  does  allow  for  preferential  weighting  and  therefore 
responds  to  recent  changes  in  trend  more  rapidly  than  the  standard  formulism. 

This  may  be  a  significant  factor  in  attempting  to  trend  a  meteorological  phenomenon 
where  the  attribute  trend  changes  from  one  value  to  another.  However,  if  the 
weighting  values  are  found  to  have  no  preferred  narrow  range  then  it  may  be  neces¬ 
sary  to  resort  to  standard  regression  methods. 

Extension  to  a  quadratic  fit  may  also  be  done.  The  standard  quadratic  relation 
is 


1'  =•  A 

tj  m 
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t+  m 
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where  this  is  usually  rewritten  as  a  linear  relation  in  terms  of  two  independent 
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or 


=Bn  + 
t+m  0 


B.  XI  4 
1  t4m 


B0  X2t 
2  t4m 


F  =  XB  . 


(17) 


F  is  now  an  N*1  vector  matrix  of  observations,  B  is  a  3*1  vector  matrix  of 
constants,  and  X  is  a  N *3  matrix  of  independent  inputs.  The  solution  is  most 
easily  obtained  through  matrix  operations  and  the  desired  constants  are  given  by 

B  =  (X'X)'1  X'Y  .  (18) 

Again  the  data  are  generally  smoothed  before  serving  as  input  to  these  relations. 
Once  again,  an  alternative  approach  which  invokes  the  use  of  weighting  constants 
and  preferential  weighting  may  be  considered.  The  Brown  quadratic  smoothing 
exponential  filter  may  be  written  as 


F.  =  At 
t4m  t 


m  4  0.5  C't  m 


(19) 


where 


At=3S;-3St4St' 


B.  =  - |(6-5cv>S:  -  (10-8o)S:'  4  (4-3a)S:"l 
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Ct=-^T-  (S'  -  2S''  4  S'") 

1  (1-ar  1  1  1 


St=aXt4(l-a)St.1 


S’’  =  aS^'  4  (l-a)S^ 


(20) 

(21) 

(22) 

(23) 

(24) 


S|"  =  aS'1  4  (l-ajs^’j  (25) 

which  as  in  the  standard  regression  method  must  be  re-evaluated  after  each  new 
data  input.  However  this  formulism  requires  far  less  computer  memory  storage 
and  computations.  Again  it  is  necessary  to  find  an  acceptable  range  of  weights 
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useful  to  the  variety  of  meteorological  environments  expected.  However,  the 
potential  difficulty  in  this  task  is  offset  by  the  realization  that  the  features  being 
tracked  and  forecast  will  generally  not  be  evolving  in  a  simple  linear  or  quadratic 
fashion.  For  a  very  short  period  of  time,  such  variation  may  appear  reasonable. 
However,  the  effects  of  using  a  finite  vertical  and  horizontal  grid  size,  coupled 
with  evolution  and  vertical  transport  effects,  and  limited  lifetime  of  the  features 
being  tracked  suggest  that  no  simple  trend  will  exist.  The  effect  of  forcing  a 
standard  linear  or  quadratic  fit  to  the  data  may  result  in  greater  forecast  error 
than  that  found  to  occur  with  these  simpler  filters  which  adapt  to  changes  more 
rapidly  and  which  may  possibly  run  with  a  generic  set  of  weights. 

One  may  consider  use  of  higher  level  forecasting  models  such  as  higher  order 
polynomial  fits,  autoregressive  moving  average  models  (ARMA),  or  Kalman  filters. 
However,  their  proper  use  requires  an  extensive  historical  data  set  for  the  feature 
being  tracked,  much  greater  than  is  expected  to  be  available  here.  Even  if  the 
data  requirements  could  be  satisfied,  these  methods  may  quickly  become  extremely 
computationally  burdensome.  Their  proper  use  may  also  require  knowledge,  or 
interpretation,  concerning  the  state  of  the  features  which  cannot  be  reasonably 
satisfied  with  the  data  at  hand.  For  example,  the  ARMA  model  requires  use  of  the 
autocorrelation  function  to  identify  characteristics  of  the  data  time  series.  The 
interpretation  of  this  autocorrelation  function  allows  one  to  choose  the  proper 
smoothing  filter  and  model.  Application  of  the  Kalman  filter  requires  substitution 
of  the  past  variance  of  the  forecast  variable  in  place  of  the  desired  current  value. 

\\  ith  the  limited  data  sets  expected  to  be  available  here,  and  the  changing  nature  of 
the  features  being  tracked,  it  is  believed  that  these  requests  cannot  be  adequately 
satisfied.  Thus,  although  high  level  forecasting  techniques  such  as  the  Kalman 
filter  may  be  tested,  it  is  strongly  felt  that  most  effort  will  be  directed  towards  use 
of  the  simpler,  more  efficient,  and  realistically  useful  methods  of  the  type  outlined 
ea  rlier. 


7.  IWMIMES  OF  METHODS 

To  bring  this  discussion  into  focus  and  provide  a  better  understanding  of  the 
proposed  solution  of  the  problem,  a  demonstration,  in  graphical  sense  will  be  pre¬ 
sented.  A  simple  evolving  data  field  will  be  studied  using  three  approaches  discussed 
earlier.  The  generalized  flow  chart  representing  the  steps  involved  in  solution  of 
tlie  problem  is  shown  in  Figure  3.  Figure  4  shows  the  three  successive  observations 
of  the  generalized  intensity  field,  which  is  portrayed  as  containing  an  advancing  line 
type  feature  from  the  west,  and  two  closed  features  in  the  eastern  portion  of  the 
viewing  area.  For  the  purposes  of  discussion  the  line  feature  is  described  by 
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simple  icon  chain  coding,  the  large  closed  contour  by  vortex  fitting,  and  the  small 
closed  feature  by  use  of  the  simple  raster  vector. 


The  three  successive  chain  codes  are  shown  in  Figure  5.  As  the  line  feature 
advances  eastward  note  that  the  chain  code  also  evolves.  Comparison  of  the  three 
codes  quickly  discloses  three  important  facts.  First,  the  main  body  of  the  code 
remains  intact,  suggesting  that  this  feature  is  basically  stationary  in  shape. 

Secondly,  the  feature  is  moving  northward  since  there  is  loss  of  code  at  the  start 
(reference  position  taken  as  most  northern  point),  and  additional  new  code  at  the 
end.  Note  that  this  is  determined  only  because  the  feature  is  believed  to  be  evolving 
slowly  in  comparison  to  our  observational  sampling  period.  Last,  rapid  structural 
change  is  localized  and  easily  detected.  The  elements  of  the  chain  code,  repre¬ 
senting  particular  locations  along  the  intensity  contour  may  be  tracked,  their  changes 
detected  and  quantified,  and  a  forecast  made.  The  dotted  lines  indicate  successive 
values  of  chain  code  for  single  elements.  The  code  along  these  dotted  lines  repre¬ 
sent  the  time  historv  of  the  individual  code  elements.  Interpretation  of  this  history 
allows  a  forecast  to  be  made  through  use  of  the  forecasting  methods  discussed  earlier. 

The  simple  raster  vector  approach  is  portrayed  in  Figure  5.  Again  a  sequence 
of  chain  codes,  with  each  code  element  representing  the  vertical  raster  line  number 
and  starting  and  ending  points  of  the  contour  of  interest  on  that  horizontal  line.is 
generated.  In  this  example  the  small  closed  contour  is  rotating  slowly  as  it  advects 
eastward,  and  splits  between  the  last  two  observations.  In  this  simple  raster 
format  the  raster  vectors  are  always  horizontally  oriented,  thus  each  raster  vector 
may  eventually  be  representing  different  contour  elements  during  the  observational 
period.  The  alternative  would  have  been  to  have  the  raster  vector  coordinate 
system  rotate  with  the  feature,  thus  allowing  the  raster  vector  to  represent  the 
same  contour  elements  for  longer  time  periods.  In  this  example  the  elements 
tracked  are  the  start  and  end  points  for  individual  horizontal  raster  lines.  The 
dotted  lines  show  the  track  sequence  history  for  individual  tracked  parameters.  It 
is  easily  seen  how  the  raster  vector  end  points  form  a  time  history  which  may  be 
modeled  to  provide  a  forecast.  The  progression  and  growth  of  the  contour  along 
the  vertical  direction  is  easily  detected  by  noting  the  changing  beginning  and  ending 
raster  line  numbers  while  the  total  separation  distance  is  essentially  maintained. 

Last,  the  splitting  of  the  contour  into  two  separate  contours  is  noted  by  the 
missing  raster  line  numbers  MB  to  20>  which  were  included  in  the  original  chain 
code  sequence.  It  is  easily  noted  that  the  primary  difference  between  this  and  the 
simpler  icon  chain  code  is  its  more  rapid  interpretation  in  spatial  terms. 

Finally,  the  vortex  fitting  approach  is  presented  in  Figure  a.  Note  that  it  is 
assumed  that  two  ellipsoids  are  used  to  fit  the  boundary  contour.  As  time  passes 
the  two  ellipsoids  are  automatically  modified  to  yield  the  best  fit  to  the  meteorological 
boundary.  The  attributes  tracked  here  are  the  ellipsoid  descriptors,  namely  the 
amplitude,  size,  center,  eccentricity,  and  orientation.  The  historv  of  each 


attribute  is  modeled  and  forecast.  The  forecast  values  determine  the  future 
ellipsoid  locations,  sizes,  and  shapes,  and  thus  ultimately  the  forecast  intensity 
contour. 
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Figure  5.  Description  of  Contour  Fields  of  Figure  4  Showing 
the  Three  Different  Attribute  Fields  and  the  Development  of 
Time  Histories  (dashed  lines)  for  Specific  Attribute  Elements 
in  (a)  Freeman  Chain  Coding,  <b>  Simple  Raster  Vector,  and 
(c)  Vortex  Fitting  Methods 


8.  CONCLUSIONS 

This  report  has  focused  on  candidate  techniques  for  forecasting  precipitation 
and  cloud  structure  at  arbitrary  locations.  A  graphical  demonstration  was  provided 
to  show  the  steps  and  applicability  of  some  proposed  methods.  It  was  determined 
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that  a  multilevel  approach  was  necessary.  Due  to  the  severe  time  constraint  of 
acquiring  new  data  and  providing  forecast  updates  every  5  to  10  min,  coupled  with 
the  reduced  error  in  tracking  and  forecasting  vertically  averaged  precipitation 
quantities  while  still  requiring  some  knowledge  of  vertical  structure,  a  simple 
three  layer  model,  representing  the  lower,  middle,  and  upper  troposphere,  is 
initially  suggested.  The  parameters  of  interest  are  the  radar  reflectivity  factor 
and  satellite  cloud  brightness.  These  parameters  require  no  additional  transforma¬ 
tion  from  their  raw  data  input  state,  and  are  readily  transformed  into  measures  of 
signal  attenuation.  The  features  to  be  tracked  and  forecast  are  the  contours  of 
these  intensity  values,  contoured  in  discrete  values.  Description  of  the  contours 
will  follow  the  process  of  isolation  by  means  of  binary  thresholding,  contour 
location  by  field  gradient  operation,  and  description  by  a  variety  of  attribute  des¬ 
cription  methods.  It  is  proposed  that  a  wide  range  of  descriptive  methods,  including 
simple  icon  and  positional  chain  coding,  up  through  complex  methods  using  poly¬ 
nomials  and  vortex  fitting  be  tested.  Simple  translation  and  rotation  of  the  whole 
contour  will  be  performed  by  cross-correlation.  Measurement  of  contour  structure 
change  was  shown  to  be  implicitly  linked  to  the  method  of  attribute  description 
through  use  of  examples.  A  range  of  forecasting  methods  were  proposed.  Although 
complex  methods,  such  as  Kalman  filtering,  will  be  tested,  focus  will  be  directed 
towards  use  of  simpler,  more  efficient  methods  such  as  linear  exponential  filtering. 
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